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ABSTRACT 

Several observational and theoretical arguments suggest that starburst galaxies may 
rival quasars as sources of metagalactic ionizing radiation at redshifts z > 3. Reion- 
ization of the intergalactic medium (IGM) at z > 5 may arise, in part, from the first 
luminous massive stars. To be an important source of radiative feedback from star 
formation, a substantial fraction (~ 1 — 10%) of the ionizing photons must escape the 
gas layers of the galaxies. Using models of smoothly distributed gas confined by dark- 
matter (DM) potentials, we estimate the fraction, (/ esc ), of Lyc flux that escapes the 
halos of spherical galaxies as a function of their mass and virialization redshift. The 
gas density profile is found by solving the equation of hydrostatic equilibrium for the 
baryonic matter in the potential well of a DM halo with the density profile found by 
Navarro et al. (1996). We then perform a parametric study to understand the depen- 
dence of {fescj on redshift, mass, baryonic fraction, star-formation efficiency (SFE), 
stellar density distribution, and OB association luminosity function. We give useful 
analytical formulas for {f e sc( z , Mjjm)) ■ Using the Press-Schechter formalism, we find 
that stellar reionization at z ~ 7 is probably dominated by small galactic sub-units, 
with M DM ~ 10 7 M and SFE ~ 10 times that in the Milky Way. This may affect the 
distribution of heavy elements throughout the intergalactic medium. 

Subject headings: Galaxies: formation - Cosmology: theory 



1. Introduction 

In all cosmological models, the temperature of the cosmic background radiation at redshift 
z ~ 1000 is sufficiently low that hydrogen ions recombine and radiation decouples from matter. 
The baryonic Jeans mass after this event is ~ 10 5 Mq, and the first luminous objects in a CDM 
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cosmology could form at redshift z ~ 30 (Peebles & Dicke 1968; Binney 1977; Rees & Ostriker 1977; 
Silk 1977; White & Rees 1978; Tegmark et al. 1996; Abel et al. 1998). If fragmentation occurs, 
triggered by H2 formation and radiative cooling (Lepp & Shull 1984; Abel et al. 1998), massive 
stars are likely to form. Depending on the slope of the initial mass function (IMF) and on the 
star- formation efficiency (SFE), the massive-star population will be a source of mechanical energy, 
heavy elements, and Lyman continuum (Lyc) photons. Such processes are known as "feedback" 
from galaxy formation, and they can have substantial impact on the intergalactic medium (IGM) 
and on future generations of stars and galaxies. 

In this paper, we focus on radiative feedback in the form of ionizing radiation from massive 
stars. In addition to its effects on diffuse gas in the halos of galaxies, Lyc radiation from the first 
stars can also provide an important source for ionizing the surrounding IGM (Madau & Shull 1996). 
Observations of the transmitted flux below the Lya emission line in high-redshift quasars (the Gunn- 
Peterson effect) imply that the diffuse IGM was already ionized at redshift z = 5. Although QSOs 
play a dominant role in photoionizing the IGM at z < 4 (Shapiro Sz Giroux 1987; Donahue & 
Shull 1987; Meiksin & Madau 1993), their dwindling numbers at z > 4 (Pei 1995; Miralda-Escude 
& Ostriker 1990; Madau 1998; Fardal et al. 1998) suggest the need for another ionization source. 
Unless a hidden population of quasars is found, radiation emitted by high-redshift massive stars 
seems necessary to reionize the universe. This hypothesis is reinforced by recent observations of 
He II absorption (Reimers et al. 1997), from the redshift evolution of the column density ratio of 
Si IV/C IV and C II/C IV in the Lya forest (Songaila & Cowie 1996; Boksenberg 1997), and form 
the thermal history of the IGM (Ricotti et al. 2000; Schaye et al. 1999) that favor a soft spectrum 
of the ionizing radiation, more typical of hot stars than of quasars (Haardt & Madau 1996; Fardal 
et al. 1998). 

A key ingredient in determining the effectiveness by which starburst galaxies photoionize the 
surrounding IGM is the parameter (f esc ), the escaping fraction of Lyc photons. This factor depends 
on the shape of the galaxy, on the gas and stellar density profiles, and on the IMF and SFE. Recent 
estimates (Madau k, Shull 1996) suggest that if 10% of the ionizing photons from hot stars escape the 
disk gas layers [(fesc) ~ 0.1], starbursts can rival QSOs as sources of the metagalactic background 
at z > 3. Haiman & Loeb (1997) made a coarse estimate of (fesc) f° r spherical objects of different 
masses at various redshift. They provided a fitting formula, log(/ esc ) = 1.92 exp[— (z — 10) 2 /i510]- 
1, valid for z > 10, and chose a value (fesc) — 1 f° r z < 10. However, these estimates are not based 
on any secure calculation or observation and we do not reproduce their results. 

Theoretical models (Dove & Shull 1994, hereafter DS94) of the radiative transfer of Lyc ra- 
diation through disk layers of a spiral galaxy like the Milky Way suggest that (f esc ) could range 
from 5-20%. The DS94 calculation considered only "burn-through": photoionized channels (H II 
regions) produced in smoothly distributed gas layers in hydrostatic equilibrium. More realistic 
models must include gas dynamics and radiative transfer through a network of radiatively ionized 
H II regions and superbubble-driven cavities ("chimneys") of hot gas (Norman & Ikeuchi 1989). By 
solving the time-dependent radiation transfer of stellar radiation through evolving superbubbles, 
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Dove et al. (2000) found that superbubbles reduce the escape fraction to (f esc ) ~ 6% because the 
shells of the expanding superbubbles can trap or attenuate the ionizing flux. By the time that 
the superbubbles of large associations blow out of the H I disk, the ionizing photon luminosity has 
dropped well below the maximum luminosity of the OB association. Observational limits on the 
transmitted flux from four low-redshift galaxies studied by the Hopkins Ultraviolet Telescope (HUT) 
are consistent with (f esc ) ranging from a few percent up to 50%. By comparing the observed num- 
ber of Lyc photons with a set of theoretical spectral energy distributions for the galaxies, Leitherer 
et al. (1995) derived, on average, (fesc) < 3%. Reexamining the same observations, Hurwitz et al. 
(1997) derived larger upper limits on the escape fraction (5.2%, 11%, 57%, and 3.2%) for the four 
galaxies. Absorption from undetected interstellar components could allow the true escape fractions 
to exceed these upper limits. 

In this paper, we attempt to provide a more realistic estimate of the production rate and 
escape fraction of Lyc from high-z galaxies. We estimate the fraction of Lyc flux escaping the halo 
of spherical galaxies as a function of their mass, integrated ionizing flux emitted by the hot star 
population, and virialization redshift of their dark-matter halos. In § 2 we describe the method 
used to estimate the escape fraction, and in § 3 we show the results of the simulations. In § 4 we 
use a Press-Schechter estimate of the distribution of virialized dark-matter halos to estimate the 
emissivity from galaxies as a function of redshift. Reionization from the first stars appears to be 
dominated by small objects with Mom < 10 7 M . Finally, in § 5 we summarize and discuss our 
results. 



2. Method 

We study the Lyc escape fraction from a spherical galaxy as a function of its virialization 
redshift, dark matter mass, baryonic fraction and total Lyc photon flux, Stot- We also consider 
the case of Stot proportional to the baryonic content of the galaxy, and we derive from simulations 
analytical expressions of (f esc ) for both cases. 

We use a Monte-Carlo method to simulate the radial distribution of OB associations in a spher- 
ical galaxy. In our fiducial model (see § 2.2), we adopt a constant stellar probability distribution 
as a function of radius, with a sharp stellar cutoff at the core radius of the dark matter halo. In 
§ 3.3 we relax this assumption, exploring the cases of a stellar density distribution that follows the 
baryonic mass distribution and the Schmidt Law. We also estimate (/ esc ) for the cases when the 
stars are located at the center of the halo and when the sharp stellar cutoff is set to a critical bary- 
onic density, ra& = 1 and 10 cm -3 . We then calculate f esc (S) for each OB association as a function 
of its Lyc photon luminosity, S (photons s _1 ), and its location in the galaxy. We calculate the 
number of OB associations with luminosities S — AS/ 2 < S < S + AS/2 integrating the luminosity 
function from the lower to the upper cutoff of the luminosity function, normalized to the total Lyc 
photon luminosity. Since we do not use a Monte-Carlo method but a simple integration, the actual 
upper cutoff of the luminosity distribution depends on the total Lyc flux. In the fiducial model 
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the luminosity function is dN a (S )/dS oc S 2 for 10 48 < S < 10 51 s _1 , but in § 3.4 we study the 
dependence of (f e sc) on the assumed luminosity function of the OB associations. Finally, we sum 
all contributions to find {f esc ) for the whole galaxy, taking care of averaging the results for several 
Monte-Carlo realizations of the same simulation. 

The method used to derive the escaping flux for a single OB association is analogous to that 
described in DS94, except that we consider a spherical galaxy instead of a disk galaxy. Thus, we 
have an additional degree of freedom, the distance R of the OB associations from the center of the 
halo. We assume that a single OB association is a point source of Lyc photons and that the gas 
is in ionization equilibrium. Since the hot stars may lie off galactic center, the shape of the H II 
region is not spherical and is determined by equating the number of photoionizations with radiative 
recombinations along each ray. Along a given ray, the Stromgren radius, r s , is defined by: 

' nij{R)r z dr , (1) 



where So is the unattenuated Lyc photon luminosity of the OB association (s _1 ), n# is the hydrogen 
number density, and a H is the hydrogen case-B recombination coefficient. We assume that the 
gas inside the H II region is isothermal at T = 10 4 K [a^ = 2.59 x 10 -13 cm 3 s -1 ]. The density 
profile of the baryonic matter is determined by solving the hydrostatic equilibrium of the gas in the 
potential well of the DM halo. We assume a spherically symmetric halo, so the density dependence 
is only on R. If the OB association lies off center at a distance D, in a polar coordinate system 
(r, 9) with the OB association centered at the origin of axes and the y-axis intersecting the center 
of the halo, we have R = (D 2 + r 2 + 2Dr sin 9) 1 ^ 2 . In Figure 1 we draw a sketch showing the galaxy 
halo with the system of reference and the notations used in the text. 

The fraction of Lyc photons escaping the halo is given by integrating over solid angle the 
fraction f(0)d0 of Lyc photons emitted between angles 9 and 9 + d9 that are not absorbed by the 
H I in the halo: 

1 fOc{S ) 

Vesc(So) = ^ J f(9)2n S m9d9, (2) 
where 9 C is the critical angle [f(9 c ) = 0] for escape (Figure 1) and 

Airrv^ f°° 

m = l-^-J Q n 2 H (R)r 2 dr. (3) 



2.1. Barionic mass distribution 

The spherical nonlinear model for gravity perturbations predicts that a just-virialized gas cloud 
has an overdensity ("collapse factor") of A c (fi m , A) = 18-7T 2 s=s 178 and is shock heated to the virial 
temperature. Thus, the baryonic gas that virializes into the DM halo of mass M^m at redshift 
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Fig. 1. — Sketch showing the galaxy halo with the system of reference and the notations used in 
the text. The OB association has a distance D from the center of the halo, and the density is a 
function of the radius R. The critical angle, 6 C , is defined for the (r, 9) polar coordinate system 
centered on the OB association. 
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z v i r has average hydrogen number density n v i r , virial radius R v i r (following convention, this is also 
called i?200i where one approximates 178 by 200), and temperature T v i r given by: 



n h h 2 \ ( 1 + z 



«„,. „ 344 pcf^V I/ 74^V / 7i±i=) , (5) 
P V 0.15 / V 10 M ©/ V 10 

O m /l 2 \ 1/3 / M DM \ 2/3 / 1 + *w \ / M 



T mr ~ 556 K — — — — - , (6) 

where Jl m , and h have the usual meanings and \x is the mean molecular weight of the gas. 

We assume a hierarchical clustering scenario for the formation of virialized halos (White & 
Frenk 1991). The results of numerical simulations show that DM halos are well fitted by a universal 
profile of mass density, p(R), valid for a wide range of halo masses. According to the results of 
Navarro et al. (1996, 1997) this profile is steeper than that of an isothermal sphere, p(R) oc i?~ 3 if 
R> R s and smoother if R < R s . Thus, we take 

£c PcO 

{R/R 8 )(l + R/R a y 

where p C Q = 3Hq/8itG is the critical density of the universe at z = 0. Integrating pdm(R), one 
finds that the overdensity is given by 

5 C = (3 x 10 3 )ft m (l + z mr f . (8) 

The characteristic radius, R S (M), of the DM is given by 



Pdm(R) = /D/D w ' , C D/D x 2 , (7) 



i 



sV ; c(M DM ) c(M DM ) \4TrA cPc0 ) ' V ; 

where c is the concentration factor, Q m is the density parameter at z = 0, A c (Q m ,A) is the 
collapse factor in a spherical nonlinear model, and z v i r is the virialization redshift of the halo. The 
concentration parameter c is related to S c by 



8, = A r 



c 3 



(10) 



_ln(l+c) - c /(l + c) 

We solve the equation of hydrostatic equilibrium of the baryonic matter numerically by assuming 
that the gravitational effect of the baryonic matter is negligible compared to the DM. The support 
of the gas in the gravitational well of the DM is dominated by bulk motions at high redshifts (Abel 
et al. 1998). We then define an effective temperature profile that is the sum of all contributions 
to the pressure (turbulent pressure, ram pressure, centrifugal force). If the gas is isothermal, the 
equations have a simple analytical solution, 

27bR s 

( R\ 

P 9 (R) = P 9 o exp(-276/2) ^1 + — J , (11) 
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where b = (2c/9r)[ln(l + c) - c/(l + c)]" 1 and where T = T vir /T g ~ 1 measures the efficiency of 
shock heating of the gas. If we assume isothermality, the resulting profile is well approximated by 
a /3-model (Makino et al. 1998), 

Pg{R) = [1 + ( R/kfW 2 ' (12) 

with [3 ~ 0.96 and R c ~ 0.22i? s . In Figure 2 we show the DM density profile and the baryonic 
matter density profile for an isothermal gas. 



2.2. Stellar density distribution and luminosity function 

We estimate the total number of escaping Lyc photons from the galaxy by integrating the 
escaping radiation of each single OB association over the luminosity function of the OB associations, 
where the position of each association inside the halo is determined by a Monte-Carlo simulation. 
We assume a power law for the luminosity function, dN a (So)/dSo oc Sq 01 for S\ < So < S*2, where 
Si and 5*2 are lower and upper limits of the observed luminosities and a = 2.0 ± 0.5, as has been 
determined for OB associations in 30 local spiral and irregular galaxies (Kennicutt et al. 1989). A 
value a = 2.12 ± 0.04 was also found for the luminosity function of young stellar clusters in the 
interacting Antennae galaxies (Whitmore et al. 1999). In our fiducial models, we assume Si = 10 48 
s _1 , 52 = 10 51 s^ 1 and a = 2; for a more exhaustive discussion of this assumption, see DS94. 

The flux emitted by all the OB associations is given by: 

Stot = fj So (^P) dSo = NobS! In (J) = JV OB (6.9 x 10 48 s" 1 ) , (13) 

where Nqb is the number of OB associations. We divide the OB associations into 10 logarithmic 
bins according to their luminosity. In the case of a = 2 and 10 logarithmic bins, the upper cutoff 
is smaller than S2 if Stot < 10 x £2. 

We assume a radial density distribution for the associations. In the fiducial models we chose a 
random distribution in radius for the distance of the associations from the center of the halo with 
a cutoff at a radius R = R s . In this region, if the gas is isothermal, the density profile is almost flat 
inside the core radius R = R c and starts to decrease near R s as R~ 3 . The average OB association 
density distribution then goes as R- 2 for R = to R = R s and is zero for R > R s . In § 3.4 we 
study the dependence of (f eS c) on the assumed luminosity function of the OB associations, and in 
§ 3.3 we study the effect on (/ esc ) of different stellar density distributions and radial cutoffs. 



3. Results 



The cosmological model we adopt is ACDM with Q, m = 0.3, VL\ = 0.7, h = 0.7 and VL^h? = 
0.019 (Buries & Tytler 1998). Unless otherwise stated, the DM halos are characterized by a 
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Fig. 2. — Profile of baryonic matter (BM, solid line) expected from the universal density profile 
of DM halo (dotted line) for an isothermal gas at virial temperature T (dashed line). The radius 
R s = i?20o/ c is related to the size of the DM core, where i?200 is the virial radius and c is the 
concentration parameter; in this plot c = 10. The density profile of the gas is well approximated 
by the isothermal /3-model with a core radius R c = 0.22R S and (3 = 0.9b (b ~ 1 is a function of the 
concentration parameter and the gas temperature). 
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Fig. 3. — Shapes of H II regions in one Monte-Carlo simulation with Nob = 100 (Stot = 6.9 x 10 50 
photons s _1 ) in a DM halo with Mdm = 10 9 M , R s ~ 1.1 kpc (outer circle), and z v \ r = 2. The 
gas core radius is R c ~ 250 pc (inner circle), and the gas effective temperature is T ~ 16,000 K. 
In this simulation we find (/ esc ) ~ 6.5%. 
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concentration parameter c = 10 and a collapse factor A c (Q m , A) ss 200. We do not consider a multi- 
phase interstellar medium, and we neglect the effect of stellar winds and supernova explosions on 
the dynamics of the ISM (Dove et al. 2000). Later in this section, we discuss these approximations. 



First, we consider a constant number of OB associations (i.e., fixed Stot), in galaxies of different 
masses. An example of a Monte-Carlo simulation is shown in Figure 3, and our results are shown 
in Figure 4. In Figure 4 (left) we show the escaping fraction as a function of the mass of the DM 
halo for a constant starburst that produces a total number of Lyc photons Stot = 3.5 x 10 51 s _1 . 
Each point is the mean of five Monte-Carlo simulations with identical parameters; the error bars 
show the variance of the mean and the three curves refer to different virialization redshifts. In 
Figure 4 (right), we show, for a fixed redshift z v i r = 10, the dependence of the escaping fraction 
on both Stot and the DM mass of the halo. The upper cutoff of the luminosity distribution is 
S U p = min(S t ot/W, S2); therefore if S tot > 10 52 s _1 it remains constant. From Figure 4 (right) we 
notice that (f esc ) remains constant if Stot > 10 52 s _1 indicating that (f esc ) depends on S U p but not 
on Stot- A good fit to the points in Figure 4 is given by a power-law function of the mass with an 
exponential cutoff at the critical mass M*, 



The lines in Figure 4 are the best fits to the points, using eq. (14) with the fitting parameters N, (3 
and M* listed in Table 1. We remind the reader that this first expression for (f esc ) is for the case 
S tot = const. 

If the DM halo is more massive than M ion ~ 6.8 x 10 9 M (l + z 1) j r )" L5 (^ m /i 2 /0.15)" 1/2 , it 
is likely that the gas falling into the potential well is collisionally ionized by the shock waves that 
virialize the gas. This effect is more likely in the outer parts of the halo. When the ionizing photons 
from an OB association reach the radius where the gas is collisionally ionized (T^,2 x 10 4 K) we 
assume that they are free to escape. Thus, the escaping fraction increases with decreasing values 
of the inner radius where the gas starts to be collisionally ionized. We show this result in Figure 5 
for the case z v { r = 10 and Stot = 3.5 x 10 51 s _1 . The curves are a parametric fit to the points, using 
cq. (14) with fitting parameters listed in Table 1. 



3.1. Case I: constant Lyc flux 




(14) 



3.2. Case II: Lyc flux proportional to the mass 



As second case, we scale the total ionizing flux, and therefore the number of OB associations 
in the galaxy, by the linear relation S to t = H*{ M g/ M g) ~ (8-75 x 10 49 s _1 )e(M fl /10 6 M ), where 
M* ~ 4 x 10 9 M is the gas mass in the Milky Way and 7* = 3.5 x 10 53 s" 1 (Bennett et al. 1994). 
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Table 1. Best-fit parameters (see eq. 14) for the curves in Figure 4 and Figure 5. 
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9.3 
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0.1 


17.0% 


0.1 


10.5 




10 6 10 7 10 B 10 9 10 10 10 6 10 7 10 8 10 9 10'° 

M DM M DM 



Fig. 4. — (left) Escaping fraction of Lyc photons as a function of the mass of the DM halo for a 
constant starburst that produces Stot = 3.5 x 10 51 s _1 . Each point is the mean of five Monte-Carlo 
simulations with identical parameters; the error bars show the variance of the mean. The four 
curves refer to virialization redshifts from z v i r = 3 to 15. (right) Same quantity for a constant 
starburst at z = 10. The three curves refer to different total number of Lyc photons: Stot = 10 53 , 
3.5 x 10 52 , and 3.5 x 10 51 s~ l . We note that (f esc ) does not depend on S to t if Sua > 10 52 s -1 . 
Some points at low mass and high total number of emitted Lyc photons shown in the figure are 
unrealistic, but we show them for illustrative reasons. 
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10 6 10 7 10 8 10 9 10 10 



M DM [M ] 

Fig. 5. — Escape fraction as a function of the dark-matter mass for a galaxy at z v i r = 10 with 
S tot = 3.5 x 10 52 s -1 . Different curves refer to situations where, outside a fixed radius cutoff R cu t, 
the gas is completely ionized by collisional ionization. 
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Fig. 6. — Results of fiducial models. Escaping fraction of Lyc as a function of the virialization 
redshift assuming a total Lyc flux, S to t = (1.14 x 10 49 s _1 )e/ s (M D M/10 6 M ), proportional to the 
mass of the galaxy, M g = fg{^lb/^l m )MDM- The error bars are the variance of the mean on several 
Monte-Carlo simulations, the dashed lines are linear fits to the mean values of (/ esc ) and the solid 
lines the analytical formula eq. (15). (top-left) Here we assume e = 4 (i.e., 4 times the Milky Way 
SFE), collapsed gas fraction f g = 0.5 and different curves are relative to M DM = 10 7 , 10 s , 10 9 M Q . 
(top-right) Assumes f g = l,e = 4 and M DM = 10 6 , 10 7 , 10 8 , 10 9 M . (bottom-left) Same as the 
top-right panel but for e = 40. (bottom-right) Same as the top-right panel but for e = 400. 
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Fig. 7. — (left) Escaping fraction of Lyc as a function of the mass of the halo for e = 4 and 
Stot oc eM;,. The solid lines are the analytical formula, eq. (15), at different virialization redshifts. 
Note the abrupt change in the shape of the escaping fraction for Stot > 10 52 s _1 when the maximum 
luminosities of OB associations reach a sharp cutoff. This behavior remain true also for different 
SFEs. (right) Escaping fraction as a function of the virialization redshift for e = 4, M = 10 7 M Q 
and gas collapsed fraction f g = 0.25,0.5,0.75,1. The dashed lines are the best fit and the solid 
lines the analytical formula eq. (15). 
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The free parameter e expresses the SFE normalized to the Milky Way value. Our adopted range 
(4 < e < 400) represents values of SFE for which appreciable Lyc escape occurs, in starbursts 
with rates much higher than that of the Milky Way. With the assumed baryonic fraction, M g = 
f g (0,b/^ m )Mr)M ~ 0.13fgMrjM, where f g ^ 1 is the product of the baryonic collapsed fraction 
and the fraction of baryons in the form of gas, S to t scales with the DM mass of the halo as 
S tot = (1.14 x 10 49 S - l )ef g (M DM /ltf M ). 

In Figure 6 we show (/ esc ) for our "fiducial models" as a function of the virialization redshift of 
the galaxies and their DM content for three different values of the parameter e with f g = l and for 
e = 4, f g = 0.5. Each point is the mean of five Monte-Carlo simulations with identical parameters; 
the error bars are the variance of the mean. In a log-linear plot (f eS c) is, with good approximation, 
a linear function of the redshift. In Figure 7 (left) we show the logarithm of the escaping fraction as 
a function of the total ionizing flux Stot (proportional to the mass Mdm) for different virialization 
redshifts at a fixed SFE (e = 4) and in Figure 7 (right) we show (fesc) for e = 4 and Mdm = 10 7 
M Q for f g = 0.25,0.5,0.75,1. We also show (solid lines) an analytical formula, 



log (/esc) 



_ 1 2 _1 

lo § (lit)'" " °- 41 ( 1 + lo S//)(^ + !) (%Sf ) 5 



lo g(l!k) 5 - 0-41(1 + log/|)( 



0.15 



,2\ 3 



1 51 s" 1 
5 2 



fgMpM 

8.8x10* Mq 



if s tot < ws 2 

if Stot > ws 2 
(15) 



that is a good fit to the points in Figures. 6 and 7 over the range of SFEs (4 < e < 400), masses, 
and redshifts shown in the figures. In the summary we write the same formula in a more compact 
form. 

The fitting formula has been derived from all the simulation results, some of them not shown 
in the paper for sake of brevity. Looking at Figure 7 (left), we note the abrupt change in the shape 
of {f esc) f° r S t ot > 10 52 s _1 , when the upper cutoff of the luminosity function stays constant. If 
we do not set an upper cutoff to the luminosity function, the dependence of (f e sc) on the mass of 
the halo is a 1/3 power-law. Instead, if we set the upper cutoff, (fesc) drops exponentially with 
the mass of the halo and (fesc) depends only on the mass and virialization redshift. In conclusion, 
small objects have bigger (f esc ) than the massive ones with the same SFE. 



3.3. Changing the Stellar Density Distribution 

The density distribution of massive stars in the halo is crucial for determining (f esc ). If all the 
stars are located at the center of the DM potential, it is trivial to find (f eS c) for a given baryonic 
density profile, because the calculation is reduced to the case of a single Stromgren sphere, 

fesc = 1 - ^ / nl( R ydr . (16) 
Jtot JO 
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(n hyo.i5)'/3(z vlr +i) (n hyo.isv^z^+i) 



Fig. 8. — Escaping fraction of Lyc as a function of redshift for different choices of the stellar density- 
distribution; the halo mass is always Mom = 10 6 M . (left) Dot-dashed lines show (f esc ) when 
all the stars are in the center of the halo for e = 4, 40. The long-dashed lines show (f esc ) when 
p* oc pi, with stellar cutoff at R(rib = 1 cm -3 ) for e = 4, 40. The dashed line has stellar cutoff at 
-R(rife = 10 cm' 3 ) and e = 4. Our fiducial models for e = 4,40 are shown for comparison (solid 
lines), (right) The dot-dashed line shows (fesc) when p* oc p& with stellar cutoff at R = R s and 
e = 4. The long-dashed lines show (f esc ) when p* oc p\ with stellar cutoff at R(nb = 1 cm -3 ) for 
e = 4,40. The other lines, shown for comparison, are the same as in the left panel. 
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Using eq. (12) for the baryonic density profile we get, 




fg n vir-R-200 > 



(17) 



and for the case Stot oc M& 



1 - 0.55/ 9 



(i+^) 3 



if (l + z)<1.22(e//,) 1 /3, 
if (l + z)>1.22(e// 9 ) 1 /3. 



esc — 







e 



(18) 



When (/ esc ) is not zero, the number of photons absorbed is such that all the gas in the halo is 
kept ionized. In Figure 8, the long-dashed lines show (f eS c) computed by using eq. (18) for two 
values of e = 4, 40. The solid lines are the fits to our fiducial models for the same efficiencies and 
Mdm = 10 6 ,10 M . It is clear that, whatever stellar density profile we assume, (fesc) will not 
be less than the value given by eq. (18) because the gas is not able to absorb more photons. In 
our calculations we do not account for the effect of overlapping H II regions. Therefore, when the 
porosity of the ionized regions approaches unity, we underestimate {f esc }. However, the collective 
effects of multiple OB associations do not become important until their H II regions overlap (i.e., 
at high porosity parameter). As shown by Dove et al. (2000), this overlap is usually accompanied 
by the production of a radiative shell, and the highest luminosity O stars have faded by the time 
this shell breaks out. Thus, we do not believe this to be a major effect. A crude correction to the 
effects of H II region overlaps would be to use eq. (18) when (1 + z) < 1.22(e// s ) 1 / 3 . 

In Figure 8 we show the effect on {f esc ) of different choices of the stellar density distribution 
and radial cutoff. A general result is that the stellar density distribution affects the dependence of 
(fesc) on the virialization redshift, while the dependence on the mass remain basically the same, 
with (fesc) decreasing as the mass increases. If the stellar density distribution follows the Schmidt 
Law (p* oc pb) and if we assume a stellar cutoff at R = R s , we find that (fesc) is exactly the same 
as in our our fiducial models, where p* oc Rr 2 and the same cutoff (solid line in Figure 8 (left)). If 
we set the stellar cutoff where the baryon density is 1 cm -3 , we find that (fesc) is approximately 
constant at high-z and increases steeply at low-z, approaching the values of our fiducial models. 
At z ~ 6, when R(nt, = 1 cm -3 ) = R s , (f esc ) is the same as in our fiducial case and approximately 
equals the constant value at high-z. This case is shown by the solid lines in Figure 8 (left) for 
e = 4,40. If we set the cutoff at R = R(rn, = 10 cm -3 ), the solid line in Figure 8 (left), the same 
reasoning holds with the only difference that the constant value of (fesc) at high-z is equal to the 
one at z ~ 15 in our fiducial case. In Figure 8 (right) we show (f esc ) assuming a stellar density 
distribution p* oc p\ and a stellar cutoff at R = R(n\> = 1 cm -3 ). The effect of a more concentrated 
distribution balances the increased cutoff radius at high-z, producing the same (f eS c) as in our 
fiducial case. 



Finally, we note that it is not unreasonable to have stars in the outer part of the halo in 
high-z objects. These halos are quite compact, and the crossing time is short compared to the 
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typical lifetime of an or B star. These stars are therefore able to move substantial distances 
from their initial location. Using eq. (5) and assuming a concentration parameter c = 10, we find 
R s = Rvir/c ~ 70 pc at z v i r = 10 and Mdm = 10 7 M . The circular velocity of a halo, defined as 
V c = [GM(R)/R} 1 ' 2 is, 

yc = (3 - 5kms ^Wj VTo^J V — To — J • (19) 

If we assume that the typical dispersion velocity of the stars is V* ~ V c , the crossing time is 
tcross = Rs/Vc — (9.6 x 10 6 yr)[(l + z 1 ,j r )/10]~ 1 ' 5 , independent of the halo mass. If the velocities 
of the gas or stars are bigger than V c , the DM halo will not be able to confine the baryonic matter 
and the galaxy will be blown away, releasing metal-enriched gas and stars into the IGM. If SNe 
explosions and stellar winds are active in this object, the stellar and gas distributions can be 
strongly influenced (Ciardi et al. 1999). In another possible scenario (Parmentier et al. 1999), the 
protogalaxy can survive to the explosion of several tens of SNe because a significant part of the 
energy released by the SNe is lost by radiative cooling. The associated blast waves trigger the 
expansion of a supershell, sweeping all the material of the cloud and the supershell is enriched 
by heavy elements. A second generation of star is born in these compressed and enriched layers 
of gas. This second generation cannot be suppressed by the Soft Ultraviolet Radiation (SUVR), 
because the cooling is provided by metals, and will have a very high {f eS c) due to the location in 
the outermost part of the halo of the OB associations. 



3.4. Changing the Luminosity Function of OB Associations 

In this section we study the dependence of (fesc) on the adopted luminosity function (LF) of 
the OB associations. We find that the choice of the LF affects the dependence of (fesc) on the 
mass of the halo. We adopt a power-law luminosity function as in § 2, but we explore the effects 
on (f esc ) of different choices of the slope a and the lower (Si) and upper (S2) cutoffs of the LF. In 
Figure 9 (left) we show (fesc) for a fixed slope a = 2 but changing the cutoffs. The solid line shows 
(fesc) when Si = 10 48 and S2 = 10 51 s^ 1 (fiducial model), the long-dashed line when Si = 10 48 
and £2 = 10 50 s -1 , and the dashed line when Si = 10 47 and S2 = 10 51 s -1 . For comparison, the 
solid line show (f esc ) for our fiducial case, (fesc) is on average lower if we chose smaller values of 
Si or S2, and it decreases more steeply with the mass when S2 is decreased. 

In Figure 9 (right) we show (fesc) fo r different slopes of the LF: a = 3 (dashed line) and 
a = 1 for S2 = 10 51 , 10 50 , 10 49 s^ 1 (solid lines). (f eS c) on average decreases if the distribution is 
steeper but the dependence on Mqm becomes shallower and vice versa. In the case a = 3 the 
relationship between (f eS c) arid Mdm is a steeper power-law. If a = 3, the total luminosity is 
the sum of many (Nob = 0.5(St o t/Si)), low luminosity associations. Instead if a = 1 the few 
(Nob = (Stot/S-z) In S2/S1) very luminous associations dominate the total luminosity. 

The simulations show that OB associations with luminosities smaller than a critical luminosity, 
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Fig. 9. — Escaping fraction as a function of the halo mass for different choices of the stellar 
luminosity function for z v { r = 6, e = 40. (left) The solid line shows (f eS c) f° r our fiducial model, the 
long-dashed line shows the effect of decreasing the upper cutoff to S2 = 10 50 s" 1 , and the dotted 
line shows the effect of decreasing the lower cutoff to Si = 10 47 s _1 . (right) The solid lines show 
{f esc) for a = 1 and (from top to bottom) S2 = 10 51 , 10 50 , 10 49 s _1 . The dashed line show {f esc ) 
for q = 3. 
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L c (Mdm), are too small to let any Lye radiation escape into the IGM unless they are far enough 
from the center. The critical luminosity increases with the mass of the halo. Therefore, in massive 
objects, only the high luminosity tail of the LF contributes substantially to (f esc ). This explains 
the decrease of (/ esc ) with the increasing mass of the halo, the dependence on the LF slope, and 
lower and upper cutoffs. When the LF is steep, only small luminosities OB associations close to 
the stellar density cutoff contribute to (/ esc ). Therefore, the dependence on the mass is shallower. 



In this section we estimate the integrated Lyc emissivity from galaxies as a function of redshift. 
The main goal is to understand the relative contribution on the emissivity from small and big objects 
and estimate the SFE (e) required to match the observed values. In our simple estimate, the number 
of photons per Mpc 3 per second emitted at redshift z is given by, 



where n(V c , z)dV c is the number density (Mpc 3 ) of DM halos as a function of circular velocity and 
redshift, M m i n (z) is the minimum DM mass for which the gas is able to cool and collapse according 
to Tegmark et al. (1996). 

To obtain n(V c , z) we use the Press-Schechter (Press &; Schechter 1974) formalism following 
White & Frenk (1991). We assume a ACDM power spectrum of density fluctuations as in Efstathiou 
et al. (1992), with Q m = 0.3, Qa = 0.7 and h = 0.7. The power spectrum normalization is 
determined from COBE measurements of the Cosmic Microwave Background (CMB) and yields 
a 8 = 1.02. 

We express the emissivity, E, in terms of the number of emitted photons per hydrogen atom 
per Hubble time. According to Miralda-Escude et al. (2000), in these units, the emissivity at 
redshifts z = (2, 3,4) should be E = (14, 8.1,4.3). These values of the emissivity have been derived 
from the observed values of the mean flux decrement at those redshifts according to Rauch et al. 
(1997). The formation of objects with Mdm < Mi on (see § 3 for the definition of Mj on ), can be 
suppressed or delayed by the SUVR produced by the first stars. The SUVR can dissociate the 
molecular hydrogen, which is the primary coolant of gas at T < 10 4 K at very low metallicity. 
Stellar feedback affects small objects by suppressing star-formation but can also increases {f esc ) by 
blowing away the gas. 

The relative contribution of small objects to the integrated (fesc) depends on the stellar lumi- 
nosity function (see § 3.4). It also depends on the number density of DM halos n(y c ,z) which is 
a function of cosmological parameters and the power spectrum of perturbations. In Figure 10 we 
show the emissivity and the mean (f esc ) (averaged over all halo masses) as a function of redshift 
(solid lines). The long-dashed lines and dot-dashed lines are the contributions to the total emis- 
sivity and mean (f esc ) from objects with Mdm < Mi on and Mdm > Mi on respectively. Triangles 



4. Simple Estimate of the Lyc Emissivity from Galaxies 




(20) 
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are the values of E from Miralda-Escude et al. (2000). In Figure 11 we show the mean (f eS c) as a 
function of redshift for different values of e and f g . 

If (fesc) does not depend on the halo mass (in Figure 10 we show the case (f eS c) = 1), the 
emissivity of the objects with Mom > Mi on and Mdm < Mi on would be equal at z ~ 13. Instead, 
because of the dependence of {fesc) on the halo mass, objects with Mdm < M{ on dominate the 
emissivity at z ^ 5 if their formation is not suppressed by feedback mechanisms. A star formation 
efficiency e ~ 8 is consistent with the observed emissivity at z ~ 4. However, if we assume that the 
reionization occurs when E ~ 1, we find z re i ~ 5. If the reionization occurs at z re i ~ 7 (Gnedin 
1999; Gnedin & Ostriker 1997) the star formation efficiency needs to be e ~ 16. Within this model, 
a constant (with redshift) star formation efficiency does not reproduce the expected emissivity 
at z = (2,3,4). In order to have a good fit to the emissivity at z = (2,3,4) and reionization 
at z re i ~ 7 the efficiency should be an increasing function of the redshift. From observations of 
interacting galaxies, Bushouse et al. (1999) have shown that the SFE increases with the amount of 
molecular gas available in the galaxy. Therefore, we expect e and f g to increase with redshift if, for 
a given baryonic mass, the gas fraction increases with redshift. Finally we remind the reader that 
our model assumes spherical galaxies and that collisional ionization of the halo gas is negligible. 
Spherical or irregular galaxies are probably numerous at redshift z ^ 2 (Dickinson 2000) but at 
smaller redshift disk galaxies are predominant. At low redshift the number of massive galaxies 
with a collisionally ionized halo increases and the effect of dust absorption also becomes more 
relevant. As a result the mean (f eS c) at z 2 should be calculated with other models (see Dove 
et al. (2000)). 



5. Summary and Conclusions 

In this paper, we have investigated the ability of the first luminous objects (Pop III stars) in 
the universe to contribute to the ionizing background that will reionize the universe between z ~ 10 
and 5. We consider a CDM cosmology, in which masses M = 10 6 — 10 7 M are able to collapse 
via H2 cooling at redshifts z ^ 30. One of the most uncertain parameters is the escape fraction, 
(fesc), of Lyc photons from the halo of the galaxies. Knowing (f eS c) is crucial for determining the 
fraction of EUV photons available to ionize the IGM. Unfortunately, the observations of very high 
redshift galaxies is still beyond the capabilities of ground-based and space-borne telescopes, and 
we can only begin to speculate about the physics of Pop III objects (Tumlinson & Shull 2000). 

However, in many theoretical models of the IGM and reionization, an estimate for the value 
of (fesc) is needed. A typical value (fesc) ~ 10 — 20% is usually adopted, based on both theoretical 
studies (Dove &: Shull 1994; Haiman & Loeb 1997) and on observations of local starburst galaxies 
(Leitherer et al. 1995; Hurwitz et al. 1997). Both estimates are valid for disk galaxies similar to the 
Milky Way at redshift z ~ 0. A recent study by Dove et al. (2000) that accounts for superbubblc 
dynamics finds (fesc) — 0.03 — 0.06. 



DRAFT: February 1, 2008 



22 




Fig. 10. — (left) Emissivity E, in photons per hydrogen atom per Hubble time as a function of 
redshift. We show (see labels in the figure) the emissivity, and the emissivity if we set (f eS c) to 
unity. The dashed lines are for halo masses M < Mi on , dot-dashed lines for masses M > Mj on 
and solid lines for all masses. The crosses are the expected values of the emissivity at z = (2, 3,4) 
according to Miralda-Escude et al. (2000). (right) The mean (f esc ) as a function of redshift (solid 
line). The dashed and dot-dashed lines have the same meaning as explained above. For both these 
figures the efficiency is e ~ 16 and f g = 1. 
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Fig. 11. — The mean (f eS c) as a function of redshift for different values of f g and e as 
from the figure labels. 
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At high redshift, the estimate of {f esc } is even more uncertain because of the lack of information 
on the SFE and IMF. For spheroidal galaxies, the stellar density distribution can also affect (fesc)- 
If the star formation is triggered by protogalaxy merging, it is plausible that many star-forming 
regions are located at the interface of the shock waves produced by the collision. If the star formation 
is not a random process, but is triggered by other star-forming regions, OB associations will tend 
to form toward the outer boundary of the galaxy, even if the first starburst happens in the center. 
In this paper, we have investigated the effect on (f esc ) of several recipes for the stellar density 
profile and luminosity function. The baryonic density profile is calculated by assuming hydrostatic 
equilibrium in the potential well of the DM halo. The gas effective temperature is assumed to be 
equal to the virial temperature, and results from both thermal and turbulent motions that support 
the gas. 

Our key results are: 

(1) The escape fraction increases with the total number of Lyc photons emitted per second, Stot, 
and decreases with increasing mass and redshift of the halo. 

(2) The stellar ionization is dominated by small galactic objects. For Stot = 3.5 x 10 51 s _1 , we 
find (f esc ) > 10% if M DM £ 5 x 10 7 , 10 7 , 2 x 10 6 , 10 6 M at z vir = 3, 5, 10, and 15 respectively. If 
we increase Stot, the DM masses listed above will increase approximatively by the same amount. 
Therefore, at redshifts z £ 10, only galaxies with DM halo masses < 10 7 M have (f esc ) > 10% 
unless the gas is collisionally ionized by shocks. 

(3) If we assume a SFE proportional to the baryonic content of the galaxy, (f esc ) decreases expo- 
nentially with the redshift and as a 1/3 power-law with the halo mass (if we set an upper cutoff 
for the luminosity function of the OB associations {f esc ) decreases exponentially for halo masses 
greater than a critical mass). The dependency of (fesc) ° n the redshift is related to the assumed 
stellar density distribution and the dependency on the halo mass is related to the OB association 
luminosity function. 

(4) We have found a simple analytical expression for (f eS c) as a function of the normalized SFE, e, 
the redshift of virialization, z v i r , and the DM halo mass, Mum'- 

log(/« c ) = log ( ) ~ °- 41 (! + l °gfg ) ( " ) (Zvir + 1) (-015" ) ' 

(21) 

where Q = S^t/lOS^. This formula is eq. (15) written in a more compact form. 

(5) A simple estimate of the emissivity as a function of redshift, using the Press-Schechter formalism, 
shows that the emissivity is dominated by small objects (Mdm ~ 10 7 M ) up to redshift 5 and 
that a SFE e ~ 8 is consistent with the observed emissivity at z = 4. A SFE e ~ 15 is needed to 
reionize the IGM at Zi on ~ 7. 

The question if these small objects with high SFE exist, has to be addressed. Their formation 
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can be suppressed by SNe explosions or by Soft Ultraviolet radiation that prevent their cooling 
via Hi. On the other hand, if they can survive to the SNe explosions from the first generation of 
stars, a second generation can born on the compressed and metal enriched gas layer produced by 
the blast waves. 

Our study is a first attempt to understand the magnitude of {fesc) from a spheroidal galaxy 
as a function of redshift. Our results are a crude approximation and can be improved by numerical 
simulations of galaxy formation. In our treatment, we do not include the effect of dust absorption, 
gas inhomogeneity, or gas dynamics. However, we believe that adding further complications to the 
model is not justified until observations tell us more about the nature of high-redshift galaxies. 
This will allow us to build a more elaborate model based on a solid observational ground. 

This work was supported by the Theoretical Astrophysics program at the University of Col- 
orado (NSF grant AST96-17073 and NASA grants NAG5-4312 and NAG5-7262). We thank Mark 
Giroux for a critical review of the manuscript and the anonymous referee for very helpful sugges- 
tions. 
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